module interp

contains

    subroutine interp_linear(nd, ni, xi, fi, x, f)
        integer, intent(in) :: nd, ni
        real, intent(in) :: xi(nd,2), fi(ni), x(nd)
        real, intent(out) :: f

        real :: t1, t2, t3, t4, t5, t6, t7

        select case(nd)
        case(1)
            f = fi(1)+(fi(2)-fi(1))/(xi(1,2)-xi(1,1))*(x(1)-xi(1,1))
        case(2)
            t1 = xi(1,2)-xi(1,1)
            t2 = xi(2,2)-xi(2,1)
            t3 = t1*t2
            t4 = x(1)-xi(1,1)
            t5 = x(2)-xi(2,1)
            t6 = xi(1,2)-x(1)
            t7 = xi(2,2)-x(2)

            f = fi(1)/t3*t6*t7+fi(2)/t3*t4*t7+fi(3)/t3*t6*t5+fi(4)/t3*t4*t5
        case(3)
            write(*,*) "3-dimension interpolation is not implemented"
            stop
        end select

        return
    end subroutine interp_linear

end module interp
